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Abstract 

We present a new formulation of the time-dependent self-interaction correction (TD- 
SIC). It is derived variationally obeying explicitly the constraints on orthonormality 
of the occupied single-particle orbitals. The thus emerging rather involved symmetry 
condition amongst the orbitals is dealt with using two separate sets of (occupied) 
single-particle wavefunctions, related by a unitary transformation. The double-set 
TDSIC scheme is well suited for numerical implementation. We present results for 
laser-excited dynamics in a ID model for a molecule and in fully fledged 3D calcu- 
lations. 
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1 Introduction 



Density Functional Theory (DFT) [Tf2f3f4"] has become over the last decades 
a widely used theoretical tool for the description and analysis of electronic 
properties in physical and chemical systems. This applies particularly to sys- 
tems with sizeable numbers of electrons [2|3] , all the more so if one is inter- 
ested in truly dynamical situations. The extension to Time-Dependent DFT 
(TDDFT) has been formally established more recently pirJfT] and it is still 
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in development, concerning both formal and practical aspects [8]. Over the 
years, TDDFT has thus become one of the few, well founded theories, allow- 
ing to describe dynamical scenarios in complex systems. This is a key issue 
for understanding dynamical microscopic mechanisms, beyond mere energetic 
considerations, e.g. the process of electron emission as it is important in con- 
nection with laser irradiation. 

A basic idea of (TD)DFT is to replace the involved correlated many-electron 
problem by an effective one-body description though the inclusion of exchange 
and correlation effects in a (as simple as possible) exchange and correlation 
functional, expressed in terms of the local density of the electrons. The sim- 
plest approximation to the exchange correlation functional is the Local Density 
Approximation (LDA), or Adiabatic Local Density Approximation (ALDA) in 
the time dependent case, which proved very useful in calculations of structure 
and low-amplitude excitations (optical response, direct one-photon processes) 
[I]. It can also be used as a first order approach in more violent dynamical 
processes involving huge energy deposits and/or large ionization as, for exam- 
ple, in the case of clusters or molecules subject to intense laser fields or to 
collision with highly charged particles [9]. 

However, LDA is plagued by a self-interaction error due to the fact that the 
direct Coulomb term and the exchange-correlation potential involve the total 
density including the particle on which the field actually acts. That Coulomb 
self- interaction is nicely canceled in a full Hartree-Fock treatment. However, 
the approximate treatment of exchange in LDA weakens this cancellation and 
a spurious self-interaction remains. As a consequence, LDA produces the 
wrong Coulomb asymptotic. The self-interaction thus spoils single-particle 
properties as, e.g., the Ionization Potential (IP) or the band gap in solids 
pUyTTj . Another critical detail where LDA fails is the polarizability in chain 
molecules [T27T3] . In dynamical situations, the self-interaction error will thus 
spoil the description of excitations involving ionization processes, especially in 
processes close to electron emission threshold. Correcting the self-interaction 
error requires a dedicated treatment known as the Self-Interaction Correction 
(SIC). Such a SIC complementing LDA static calculations was proposed in 
[T^fTo] . It has been used since then at various levels of refinement for struc- 
ture calculations in atomic, molecular, cluster and solid state physics, see e.g. 
[T6fl7|ll8|rT9] . The original SIC scheme, however, leads to an orbital depen- 
dent mean field which causes several formal and technical difficulties. This 
aspect can be circumvented by treating SIC with optimized effective poten- 
tials (OEP), see [20] for a recent review. The resulting formalism is quite in- 
volved and usually treated with involving further approximations, as e.g. the 
Krieger-Li-Iafrate (KLI) approach [2TH22"] . These, however, can perturb some 
crucial physical features of SIC, particularly the trend to produce localized 
single-particle states [20] . 
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The applications of SIC in time- dependent situations have, up to now, mostly 
been performed in the above mentioned approximate manners, e.g., the lin- 
earized treatment of [23], the use of averaged-density SIC [21], or the various 
versions of time dependent OEP-KLI [2Bl26f27j . The latter TDOEP-KLI, how- 
ever, also suffers from inconsistencies. It leads, in particular, to the violation 
of zero force theorem and energy conservation [28J . There thus remains plenty 
of space for elaborating adequate versions of SIC, and even more so of TDSIC. 
The natural starting point, and benchmark for later approximations is a full 
TDSIC scheme. The aim of this paper is to present an exact thorough varia- 
tional formulation of fully fledged TDSIC. We shall also propose a manageable 
propagation scheme which allows to obey all key conditions, namely the zero- 
force theorem, conservation of energy and conservation of orthonormality of 
the occupied single-particle orbitals. Our approach relies on a simple account 
of basic constraints and on the use of an important degree of freedom, namely 
the freedom of unitary transforms among occupied orbitals. 

The paper is organized as follows. We first remind basic SIC equations in 
static case, introducing already the unitary transform degree of freedom. We 
then export the formalism in the time domain and discuss the properties of 
TDSIC. We finally show practical examples of applications in simple molecules 
and clusters, in particular in the case of irradiation processes. 



2 Basic notations 

Before attacking the question of the self-interaction correction, we want to 
introduce briefly generic notations and take the example of widely used Lo- 
cal Density Approximation (LDA) [29], [30] which serves as a basis for our 
further considerations. We shall work in the Kohn-Sham scheme of DFT 
|31j . The Kohn-Sham state is composed of a set single-particle wavefunctions 
{ip a , a = 1, . . . , N} where N is the number of electrons of the system. These 
single-particle states have to be orthonormalized. This requirement will play a 
role later on. Both static and dynamical DFT schemes then amount to write ef- 
fective one-body Schrodinger-like equations for the {ip a) a = 1, . . . , N}, called 
Kohn-Sham (KS) equations. In the KS scheme, the total electronic energy of 
the system E can be split into four terms : 

E = E kin + E cxt + En + E xc . (1) 

The kinetic component E^ in is computed assuming non interacting ip a \ the 
direct Coulomb interaction En is computed computed classically (Hartree ap- 
proximation [32]); the effect of the external potential (E cxt including in par- 
ticular the ionic potential and possibly external fields such as that delivered 
by a laser) is computed exactly and finally the exchange Coulomb and the 
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electronic correlations are packed into the exchange correlation energy E xc for 
which one has to construct approximations. All terms are functionals of the 
total electronic density p(r). 

The simplest and most widely used approximation for E xc is the Local Density 
Approximation (LDA) in which one performs a local Fermi gas approxima- 
tion for evaluating energies. The LDA serves as a starting point for many more 
involved approximations, in particular the SIC approximation we discuss in 
this paper. We thus assume that E xc is computed in the LDA approximation. 
For the sake of simplicity in the notations, we shall pack together the (exact) 
Hartree and exchange correlation terms and note -Elda the corresponding en- 
ergy at LDA approximation : -Elda = En + E xc . One can then derive the 
KS equations (stationary or time dependent) by standard variational tech- 
niques which leads to single electron KS equations with LDA single electron 
Hamiltonian 



K 2 A 

h~LDA = - -z h U cxt + U LDA [p] , (2) 

2m 



^lda[^] 



SE LD a 



Sp 



Uh[q] + U xc [q] . (3) 



P=Q 



As outlined in the introduction, the LDA approximation suffers from the self- 
interaction error. For example, its direct Coulomb part Un[p] is a functional 
of the total density which is computed by summing over all occupied single 
electron densities (p = J2 a Pa)- Then a particle a will feel its own Coulomb 
repulsion and this spurious self-interaction is not properly removed by the 
exchange term in LDA. This thus calls for a SIC treatment. 



3 Stationary SIC 

3. 1 SIC functional and Hamiltonian 



The starting point is the SIC energy functional (following notations of section 



N 

E S ic = E kin + E cxt + E ion + E LDA [p] - E LDA Uf3 1 2 ] (4) 

where the total density is p = J2a Pa with p a = \i/j a \ 2 - Note that all summations 
run over occupied states only. The corresponding one-body Hamiltonian is 
obtained from variation of E$ic with respect to ^* as 
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— — = h a ip a , (5a) 
6iff* 

K = h LDA - U a , (5b) 

U a = Ujj,M a \ % \ . (5c) 

The emerging one-body Hamiltonian h a depends on the state ip a on which it 

acts through the SIC term U a . Thus it is not invariant under unitary trans- 
formations within the sub-space of occupied orbitals. 



3.2 The stationary SIC equations 



3.2.1 Variational derivation 

The static SIC equations are derived by minimization of the SIC energy (j4j) to- 
gether with the condition that the single-particle orbitals are orthonormalized. 
This amounts to the variational equation 



= cL* 



£sic-E(^)A 



a,P 



(6) 



where X a p is a matrix of Lagrangian multipliers, which is non-diagonal in 
general. As worked out in appendix |A] \ a p is a hermitian matrix. Evaluation 
of the variation yields the stationary equations 



X/3a = ('ipl3\h a \lp a ) . 

We consider the hermitian conjugate equation 



(7a) 
(7b) 



(where we have exploited hermiticity of \pa), project Eq. ([Taj) with (^l, its 
conjugate with \ip a ), and take the difference of these two equations. This yields 
= (ijjp\hp — h a \ip a ) . The only state-dependence in h a stems from U a . Thus 
we remain with the condition 



0= (i>p\Up - U a \i/; a ) 



(8) 



We call it the symmetry condition. It plays a crucial role in all SIC consider- 
ations. It was first introduced in a particular case by Pederson et al. [16] and 
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since then addressed by several authors [T7|19ll29j . The above derivation indi- 
cates clearly the relation between symmetry condition Eq. (jHJ) and orthonor- 
mality constraint. We shall discuss this condition further at several places. 
Note that it is trivially fulfilled in case of state-independent Hamiltonians for 
which h a = hp = h, as it is the case in LDA or Hartree-Fock. 

A word is in order about the SIC Hamiltonian h a . It depends on the state on 
which it acts. Thus one has to be extremely careful with everything one knows 
from Quantum Mechanics and Hilbert space. The h a is not a linear operator 
which is obvious from the fact that the operation h a [|^> a )c a + \ij}p)cp\ is not 
defined at all. We will also see more clearly in the next section that the SIC 
Hamiltonian is not hermitian. 

3.2.2 State-independent notation 

For formal manipulations, it may be simpler to recast the state dependent SIC 
Hamiltonian h a into a compact form as 

^SIG = ^LDA-5^^«I^«)(^«I ' ( 9 ) 
a 

The part sensitive to single-particle states has been expressed in terms of 
projectors \ip a )(ip a \ such that the SIC Hamiltonian (Q is not explicitly state- 
dependent, but hsic\ipa) remains equivalent to h a \ip a ). The form ([9]) is advan- 
tageous for formal considerations. The SIC equations become now equivalently 



^SIC#a) =J2 I^V* > ( 10a ) 

p 

A/3a= {i>p\hsic\^a) ■ (10b) 

The symmetry condition is derived as above. We build the hermitian conjugate 
equation, take the difference of the two equations, and exploit the fact that 
the Lagrangian matrix is hermitian. This yields 

= {^Mic ~ kicM (10c) 

which is equivalent to the symmetry condition in the form (JSj) when dropping 
the state independent part /zlda and evaluating the projectors. 

It is interesting to note that the form (Q shows clearly the possible non- 
hermiticity of the SIC Hamiltonian. This makes the symmetry condition (I lOcj) . 
or equivalently (JS}, a non-trivial and crucial part of the SIC equations. Su- 
perficially, it makes the impression of a condition ensuring hermiticity of the 
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SIC Hamiltonian fasic- But one has to keep in mind that condition (1 10c[) is re- 
stricted to the space of occupied orbitals. Thus the symmetry condition forces 
restoration of hermiticity only in the sub-space of occupied states. 

3.2.3 Projector notation 

The SIC equation can be recast in a particularly compact form when intro- 
ducing the projection operator onto the unoccupied space 

n± = i-EliWWv»l • (ii) 

& 

This allows to reformulate the SIC equations (ITaj) or fllOal) as 

n±^a) = o or nAicl^a) = o (12) 

showing that these equations serve to establish a decoupling of occupied and 
unoccupied space which is a general feature of any mean-field equation. The 
new key feature of the SIC equations is the additional symmetry condition, 
Eqs. (IH]) or (II Pel) , which comes into play because the SIC energy (jlj) is not 
unitary invariant such that there is a unique optimum for the occupied states. 

3.2.4 Single-particle energies 

The symmetry condition minimizes the SIC energy and does that by producing 
more or less localized states which maximize the Coulomb SIC of each state 
(see the later discussion). This produces in general non-diagonal Lagrangian 
matrices \ a p from which single-particle energies cannot immediately be read 
off. However, the necessary information is contained in that matrix. The single- 
particle energies can be defined as the eigenvalues e» of \ a p obtained from the 
secular equation in occupied space 

^apVpi = ZiVai (13) 



where v a i are coefficients of the appropriate unitary transformation. 
3.3 Double- set formulation of SIC 

Thus far, the formulation of SIC for stationary states is complete and man- 
ageable. The computation of single-particle energies motivates an alternative 
formulation which deals with two different, but related, sets of occupied single- 
particle states. We will thus discuss in this section a double-set formulation 
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of stationary SIC. It is an interesting, but not compulsory, alternative for the 
static case. But a double-set technique becomes almost inevitable for TDSIC. 
The present (static) section serves, so to say, as a preparation. 



3.3.1 Two sets of occupied states 

The computation of single-particle energies, as outlined in section I3.2.4[ leads 
naturally to a second set of single-particle states {<fii} connected to the original 
set by a unitary transformation within occupied space 

N 

¥i=Yl ^^Vai , J2 V ai V »j = $ij ■ (14) 
o=l a 

The set {<fii} is associated to the single-particle energies £j and so diagonal 
in energy space. We call it diagonalizing set . The set {ip Q } optimizes the SIC 
potentials and does that by some localization. We call it the localizing set . 
The diagonalizing set is compact in energy space at the price of larger spatial 
spreading and the localizing set minimizes spatial extension while enhancing 
energy variance. Both sets have their value. A proper combination of them 
will become particularly important in the dynamical case, see section HJ 



3.3.2 Double-set SIC equations 

The first SIC equation fllOap becomes particularly simple in terms of the di- 
agonalizing set. It reads now 

h<sic\Wi) = ei\<pi) . (15a) 

That equations provides the decoupling from unoccupied space as shown in 
section 13.2.31 The symmetry condition cares for determining the localizing set 
within occupied space which now shrinks to a condition for the transformation 
coefficients v ai . We emphasize that by rewriting 

Via <— > = tyf,\Up - Uatya) • (15b) 

That localizing set is needed to compute the SIC potentials U a and with it, 
ft-sic- 

The double-set equations (fT5j) can be used for an alternative solution scheme. 
However, there is no gain in efficiency as compared to the previous scheme, 
i.e. solving first the SIC equations (0) with a single set {ipa} and afterwards 
diagonalizing the Lagrangian matrix \p a to obtain the single-particle energies. 
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3.3.3 Variational derivation 



It is instructive to derive double-set SIC directly from the stationary varia- 
tional principle. To that end, we consider the diagonalizing set {y^} and the 
transformation coefficients v a i as variational degrees of freedom. The SIC func- 
tional is to be minimized with boundary conditions of orthonormality of the 
tfi and v a i. This means to minimize the functional 

F[<fii, V^] = E S Ic[<fi, Vai] ~ Y,(Vk\¥j)9jk ~ ^ (l2 V ai V (3i) A /3a ■ (16) 

k,j otfi i 

It is interesting to compare the extension with LDA. In that case, one deals 
with an energy functional which is invariant under unitary transformations 
amongst occupied states. That allowed to perform always a unitary trans- 
formation such that ^kjifklfjWjk — ► ^j(fj\fj) £ j fr° m which one obtains 
immediately the energy-diagonal LDA equations by variation. The SIC func- 
tional (j3]) is not unitary invariant which, in turn, led to the notoriously non- 
diagonal Lagrangian matrix. The functional (fTBT) formulated in terms of the 
double-set can now be considered again as being unitary invariant with respect 
to the ipi because any rotation within the {<fi} can be compensated by proper 
counter-rotation of the Vi a . Thus we can always perform a transformation to 
the simpler functional 

j a/3 i 

First, we perform variation with respect to the <p*. The key piece is 



= (r\h$ic\<Pi) ■ 

Thus we obtain from d^F = the first SIC equation ( 11 5 aft . In a second step, 
we perform variation with respect to the transformation coefficients v ai . We 
exploit 

SEsic f-,3 Si) a 5E S ic . .f . v i \t \ i \ 



Sv ai J 5v ai 5ip a 
and obtain from variation J2p Vpifyplhalipa) = J2p J2i v pi A pa and subsequently 
(ipp\h a \x/) a ) = Ap a . Similar as in sections 13. 2. Il and l3. 2.21 we build the hermitian 
conjugate and take the difference. This yields then the symmetry condition 
(HIS). 

Thus the direct variational derivation recovered nicely the double-set formula- 
tion of stationary SIC. It is important to remark that the symmetry condition 
(115bp results from minimization of the SIC energy with orthonormality con- 
straint for fixed </?j. 



9 



3.3.4 The existence of a solution to the symmetry condition 



The symmetry condition Eq. (115bj) is as such a highly non-linear equation. 
One may wonder whether a solution exists in general. We have seen in the 
above section that the symmetry condition simply emerges from minimizing 
the total energy in the reduced space of occupied single-particle orbitals. There 
necessarily exists an energy minimum in the restricted space and thus that 
there always exists a solution to the symmetry condition. This is a crucial 
feature because the symmetry condition is always present in any formulation 
of SIC, static and time-dependent. 



4 Time Dependent SIC (TDSIC) 



Now that a proper SIC formulation has been given in the static case, we can 
consider the dynamical case along the same line. The diagonal formulation 
will become crucial. 



4-1 Derivation of TDSIC 



The TDSIC equations are obtained from the principle of stationary action 
using the SIC energy functional fll]) 



o = ss , s= [Ut'(E slc -J2MmM-^(Mi>j)^) , (17) 

J to \ a p a / 

explicitly including the orthonormality constraint with Lagrange multipliers 
A 7j g as in the static case. Note that the matrix of Lagrangian multipliers is 
hermitian as shown in appendix E Variation with respect to yields the 
TDSIC equation for the propagation of single-particle orbitals as 



(/isic - iftdt) \ij) a ) = ^ IVV?) V* , (18a) 
P 

A/3a = (ipp\h a - ihdt\ip a ) - (18b) 

The relation f)18b[) for the Lagrangian multipliers becomes non-trivial by the 
fact that it is hermitian, i.e. 

A/3 Q = A*^ . (18c) 

This can be exploited by the same steps as performed in the static sections 
13.2.11 and 13.2.21 We build the hermitian conjugate of Eq. (118bl) . insert Eq. 
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(I18cp . and take the difference. This yields once again the symmetry condition 

= (ii>f,\Uf,-U a \il> a ) , (18d) 
now for TDSIC and to be fulfilled at each instant of time. 



4-2 Solution of TDSIC with a double-set of orbitals 

The TDSIC equations (|18p are very involved and it is extremely hard to deduce 
a transparent numerical stepping scheme from them. Time evolution is related 
to energies and we have seen in static SIC that single-particle energies and 
SIC potentials are taking different cuts through the single-particle Hilbert 
space. The concept of single-particle energies led us naturally to a double- 
set strategy, see section 13.3.21 That strategy becomes extremely helpful in 
developing a solution scheme for TDSIC. 

We disentangle the involved equations of motion (1181) by distinguishing the SIC 
localizing set {ip a (t)} from a propagating set {<pi(i)}. The both are connected 
by a unitary transformation amongst occupied states 

N 

\<Pi{t)) = E \M*)) "fa® , E<*(*)<;(*) = %- (is) 

That is the time-dependent generalization of the transformation (fl4|) . The 
transformation coefficients depend also on time and the transformation is per- 
formed at each instant of time. 

We now choose the propagating set <fi such that it diagonalizes the Lagrangian 
matrix \ a p. Thus we obtain (hsic ~ ifidt) Wi) = -^mIv 9 *)- The yields an 
irrelevant phase and can be ignored. There remains 

(/isic - ihd t ) \<Pi) = , (20) 

The {y?i} can then be propagated in standard manner as: 

\<Pi{t)) = exp j~ £ dt'h sl c(t')}\<Pi(t )) • (21a) 

This procedure implies that the symmetry condition (118d[) is fulfilled at each 
instant of time. To that end, we exploit the freedom of choice of v a i. Similar 
as in static SIC, we know that we can always determine the v a i such that 

V ai (t) <— > 0= tyfilUp - U a \l/> a ) . (21b) 

The interlaced stepping of Eqs. (I2TT) provides a manageable solution scheme. 
Further formal properties will be discussed in sections 14.41 and 14.51 The prac- 
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tical applicability and stability of the scheme is proven by the many results 
presented in section |5j 

Note that the above propagator in Eq. (12 lap is not strictly unitary because 
hsic is not hermitian. But the hermiticity within occupied space, Eq. (118c|) . 
guarantees that the propagation (12 lap preserves orthonormality within occu- 
pied space, i.e. (<Pi(t)\<Pj(t)) = That suffices for our purposes. 

4-3 Projector notation 

An instructive alternative formulation can be given using the operator of pro- 
jection onto the unoccupied space, the II j_ as defined in eq. (jlip . This allows 
to recast the TDSIC equation into the particularly compact form 

n±(m - hsicMa) = . 

It defines the part of the change of the wavefunctions evolving into the space 
orthogonal to the already occupied states. The evolution inside the occupied 
states is again prescribed by the symmetry condition (j21bj) . 

4-4 Direct variational formulation of double- set TDSIC 

In section 14.11 we deduced TDSIC from variation of the action (ITTj) with 
respect to the single-particle states {ip a }, in a one-set strategy, while imposing 
their orthonormality. The double-set TDSIC was introduced in section I4~2l as 
a means to solve the TDSIC equations. In this section, we are going to derive 
double-set TDSIC directly from the time-dependent variational principle. This 
makes the derivation of TDSIC especially straightforward and it will add new 
aspects to the scheme. 

Starting point is again the action (ITTj) . but now formulated in terms of the 
two sets of orbitals, the SIC orbitals {ipa} and the propagating orbitals {<fi}, 
related by a unitary transformation ()19p . Variation with imposing orthonor- 
mality of the {ipi} and {v ia } reads 

S[ipi,V ia \ = I 6£ (Esic[^a\ - X^I^N) _ 

Jt ° ^ a k,l 

-E(E<^)V) • (22) 

a/3 i J 

As proven in appendix [Xj, the Lagrangian matrices 8jk and A^q, are hermitian. 
Note that the transformation (I19p leaves the time- derivative term invariant 
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and we have chosen to express it in terms of the propagating set which is here 
the natural choice. The SIC energy (jlj) is not unitary invariant and needs to 
be expressed in terms of the SIC set {ip a }- We ought to remind, however, that 
the {4> a } are given through the {<pi} via Eq. <HM such that we consider the 
action as a functional of the {<£>?,} and {v a i}. 

First, we perform variation with respect to the coefficients of the unitary 
transformation. We note that, among the first three terms, only E^ic depends 
on them and thus S v S = leads to 

8v ai ( ^SIC V L^i(3 A l3a ) = 

V a/3 i J 

and subsequently to 

/ A (g§-?***)-° ■ 

where h^ic is given in Eq. (Q . That can be rewritten in the more familiar form 
as (<Pi\h a \i/} a ) = J2p v pi-Ap a and finally be transformed to 

(i/)p\h a \i/) a ) = Ap a . 

Considering the complex conjugate of that equation and exploiting hermiticity 
of Ap a , we find (i()p\hp — h a \ip a ) = and from this the crucial symmetry 
condition 

v ai < — > {^pWp - U a \i> a ) = , (23a) 

which is here to be understood as a condition determining the coefficients Vi a 
for given set {tpi}. 

In a second step, we perform variation with respect to the propagating orbitals. 
Evaluating the variational equation S v * S = yields 



(hsic - ihdt)\<pi) = Mfyi > °ji = (vAic - i^lVt) 

3 

The Lagrangian matrix is again hermitian, i.e. 8^ = 9*^ which, in turn, implies 
the "weak" hermiticity condition that /isic is hermitian in the sub-space of 
occupied states. Thus this Hamiltonian can be diagonalized and it is sufficient 
to solve (/isic — ihdt)\<Pi) = Vil^i)- The Floquet index r]i produces a global 
phase factor which is irrelevant for our purposes and can be dropped. Thus 
we remain with the time-dependent mean-field equation for the propagating 
states 

(h slc -m)\<Pi) = ■ (23b) 
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The propagation (123bj) together with the symmetry condition (I23al) constitute 
the complete set of dynamical equations for TDSIC with double set. It is 
satisfying to see that both equations can be derived out of one variational 
principle. The difference to TDSIC with one set as derived in section 14.11 is 
that we now allow for two independent sets of orbitals connected by a unitary 
transformation. The variational scheme exploits that additional freedom to 
deliver correctly double-set TDSIC as the optimal scheme when dealing with 
two sets. 

We see also that the symmetry condition emerges from variation of -Esic much 
similar as in the static case. The reasoning of section 13.3.41 proving the exis- 
tence of a solution for the symmetry condition does also apply here. 



4-5 Conservation laws 

The TDSIC equations yield energy conservation as long as the external field 
remain independent on time. And they also fulfill the zero-force theorem 
[2"5|33|34j . The reasoning is simple. The LDA functional and its SIC exten- 
sion are invariant under time- and space-translations. The same holds for the 
orthonormality constraints. The equations of motion are derived variationally 
without further restrictions and approximations. This yields energy conser- 
vation from time translational invariance. The space-translational invariance 
would yield momentum-conservation if the electrons were alone. That feature 
is broken by external fields. But what remains is the fact that the electrons 
cannot exert a force on themselves which is the content of the zero-force the- 
orem. More explicit proofs are given in appendix [D] 



4-6 Static limit 

It is also interesting to consider the stationary limit of TDSIC. To that end, 
we use TDSIC in double-set set formulation. We identify 

^(r,t) = ^(r,0)e^'/ a . 

Inserting that into the TDSIC equation (T2"U1) . one immediately recovers the 
static SIC equation f)15ap for the diagonalizing states. The single-particle wave- 
functions change in time only by a phase factor. The sub-space of occupied 
states thus remains constant in time and the symmetry condition always min- 
imizes the same sub-space. And that yields the time-independent localized 
states ij; a of the static problem. 
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5 Test cases and numerical realization 



5.1 Models 



We want to test the above discussed TDSIC formalism in truly dynamical 
situations. The ultimate goal of these studies is to describe dynamically the 
irradiation of various molecules, in particular organic ones. A correct modeling 
of the IP is then crucial and TDSIC becomes compulsory, especially close 
to emission threshold. We have developed since long fully fledged coupled 
electronic dynamics for clusters [35] and used extensively at TDLDA level 
and in simplified SIC schemes [21]. We thus have at hand a powerful tool for 
analyzing irradiation dynamics. Still, the proposed TDSIC formalism is quite 
involved. In order to have a more flexible testing tool, we have thus developed 
a one-dimensional (ID) model mocking up typical atoms and simple (linear !) 
molecules. This will serve as a starter to study the properties of TDSIC. And 
we will finally complement these schematic results by realistic ones with the 
full 3D approach. 



5.1.1 The simplified ID model 

To test numerically the SIC scheme, we take up the test case of [36] consisting 
in a one- dimensional model for a molecule. Spin is not taken into account 
explicitly and all electrons are assumed to have the same spin such that they all 
explore the full exchange effects. Apart from its simplicity and computational 
cost, the ID test case has also the advantage to be a much more sensitive test 
of orthonormalization than in 3D calculations. 

For the electron-electron interaction, we use a smoothed Coulomb potential 
(in Hartree units) 



Starting with this "elementary" interaction, we construct the corresponding 
LDA energy functional for exchange only. Working at the level of exchange 
only allows to have fully fledged time-dependent HF (TDHF) calculations as 
a benchmark to which TDSIC calculations can be compared. The detailed 
calculation of the LDA energy is presented in appendix O The resulting LDA 
exchange potential (7 is the possible degeneracy number, equal to 1 in the 
next results) reads: 




1 



(24) 





2np(x) 



(25) 
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For the electron-ion interaction, we also use a smoothed Coulomb potential 
(in Hartree units) of the form 



Nz Ml - z) 

U ion {x) = / - 1 ; , 26 

(x - R/2) 2 + b J(x + R/2f + b 



where R is the inter-ionic distance which allows the possibility to compute bi- 
atomic molecules, N is the number of electrons and z G [0, 1] is an asymmetry 
parameter. If z = 0.5 and R=0, we recover the atomic case. If z=0.5 and R ^ 
0, we obtain the bi-atomic symmetric molecular case. If z ^ 0.5 and R ^ 0, 
we get the asymmetric molecular case, which provides a more critical probe 
of the role of orthonormalization. 



The a and b parameter for a two-electron system are scaled to approach the 
experimental H 2 bond length (1.4 ao) and ionization potential (0.57 Ha). With 
a=0.8 and 6=0.5, we obtain a H 2 bond length of 1.6 a and an ionization 
potential of 0.5 Ha. Depending on the case, we use slightly varied values which 
will be indicated. 



5.1.2 Full 3D model 



The 3D calculations follow the standard techniques as we use it since long 
[91135] . The electron wavefunctions and spatial fields are represented on a 
Cartesian grid in three-dimensional coordinate space. The spatial derivatives 
are evaluated via fast Fourier transformation. The ground state configurations 
were found by adapting the accelerated gradient iteration for the electronic 
wavefunctions [37] to SIC. Propagation is done by the time-splitting method 
for the electronic wavefunctions [38] augmented by updates of the symmetry 
condition as explained below. For the energy functional, we employ the widely 
used functional of |39|. 



5.2 Solution scheme for stationary SIC 



In straightforward generalization of the damped gradient step [37], we solve 
the static Eqs. (J7J iteratively as 



^stcp 



T + E, 



damp 



P 



(,ipp\h a \ip a ) + (ip a \hp\ipp)* 



(27a) 
(27b) 
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where O stands for orthonormalization. The step (I27aj) provides simultane- 
ously a complete solution including a matching of the symmetry condition, 
because the latter is implicitly taken into account in the symmetrization of 
the A matrix. 

It turns out that the step (I27al) converges, however very slowly, and that the 
symmetry condition causes the delay. We speed up the iterations improving the 
symmetry condition explicitly in each step. This is done by a unitary transfor- 
mation within occupied states. The coefficients of that unitary transformation 
are also determined by a gradient iteration as 



D 



(old) 

47 



77 A 7 } 



d u * 



E S IC ~ U *ja U jt3^f3a 



-Oil £^7) -J2 u i/3 X ^ 



(28a) 
(28b) 



where the "driving force" is obtained by variation of the SIC energy 
with respect to u* . That interlaced combination of damped gradient step and 
symmetry condition converges acceptably fast. 

Depending on the initial conditions, it may take a while until the interlaced 
iteration has found its path to the properly localized wavef unctions. A further 
substantial acceleration can be achieved by performing in the initial phase once 
in a while a localization transformation. There are several localization criteria 
at choice P0|^T] . We found very efficient improvements with simply minimizing 



the sum of the spatial variances of the single-electron states, defined as 



M = E [(^|r 2 ivg - (^*|i# 



(29) 



5.3 Propagating TDSIC numerically 



In the time- dependent case, the only manageable way to propagate TDSIC is 
to use the double-set strategy. For one time step St_, the propagation proceeds 

as follows. We first evaluate ifi (t + — ] = exp 



iSt f . . 
-2h hsic{t) 



fi(t), where 



hsic{t) is obtained by the chain : 



<Pi\t) ► »■ h S ic{t). 



(30) 



Gradient iteration, similar to Eq. (128b|) . is used to solve Eq. (jl8d.fl for the 
Vip(t), from which one deduces the if) a (t) and hsic(t). The <fi(t + St/2) thus 
obtained are used to compose hsic(t + St/2) similarly using the chain (13DI) . We 
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finally compute 



<Pi(t + St) = exp 



i6t 



h 



SIC 



t + 



5t' 



<Pi(t). 



(31) 



After all, the scheme as explained here works reliably and robust in ID as well 
as in 3D. We checked conservation of energy and orthonormality and found it 
fully satisfying in all cases. The symmetry condition is, of course, fulfilled all 
along the time evolution by construction. 

A proper computation of ionization requires to prevent reflection of electrons 
which have been ejected from the molecule and are now impinging on the 
bounds of the numerical grid. We do that by employing boundary conditions. 
To that end, an absorbing zone of a few grid points is defined. For each time 
step, the part of electronic wavefunctions which has penetrated into that zone 
is eliminated by a mask function, for details see [9|35] . 



6 Results and discussion 



6.1 Stationary state 



0.2 

o 

■0.2 
■0.4 
■0.6 
■0.8 
-1 
■1.2 
-1.4 



Ionic + LDA potential 



-2 2 4 6 
x (Bohr) 



Fig. 1. The mean- field potential from ionic background and LDA-part (without 
SIC-parts) obtained in the ID model with 2 electrons. The model parameters are 
a = 0.8 ao, b = 0.5 ao, R = 1.5 ao, and z = 0.4. 

It is well known that static SIC has a tendency to localize spatially the orbitals 
[T6] . We analyze this fact on the example of the stationary solution for a system 
of two electrons in the ID model developed in section 15.1.11 Note that z ^ 
0.5, i.e. we enforce a slight asymmetry. The resulting ionic + LDA potential 
obtained for a SIC solution is plotted in Fig. CD Note that the SIC potential 
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cannot be plotted easily because of its state-dependence. The resulting single- 
particle energies e, are -0.88 / -0.32 for LDA, -1.18 / -0.60 for SIC, and 
— 1.24 / —0.55 for the HF benchmark. Note that the results of SIC comes close 
to HF as it should be. 




-6 -4 -2 2 4 6 -6 -4 -2 2 4 6 

x (Bohr) x (Bohr) 




-6 -4 -2 2 4 6 
x (Bohr) 

Fig. 2. Static SIC densities obtained in the ID model with 2 electrons. Top left: 
Densities from diagonal wavefunctions, \(fi\ 2 . Top right: Densities from localized 
wavefunctions, (V'al 2 ) obtained after unitary transformation. Bottom: Comparison 
of the total density calculated from both sets of wavefunctions. 

Fig. [2] compares the single-particle densities of the localized wavefunctions 
| i/j a | 2 with those of the energy-diagonal wavefunctions \tpi\ 2 . It is obvious that 
the localized densities are much better concentrated in space than the diagonal 
ones. The total density (lower panel), of course, remains the same for both 
cases because the two sets are linked by a unitary transformation. 

We are now considering dynamical evolution from the given static state, where 
the excitation is initialized by a very short laser pulse, simulated as an instan- 
taneous boost [12]. The question is to what extent the localization may be 
washed out through the excitation. In order to follow an evolution, one needs 
to characterize localization by one number and we do that by the spatial 
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Fig. 3. Time evolution of single-particle variances (j32f) obtained in the ID model 
with 2 electrons and asymmetrical ionic background. Shown is the relative value 
(Atp — Aip)/Ail)). The time averaged result is indicated by a straight line. 



variance of a single particle states 



x 2 \ipi 



^(</?i|x 2 |<^) 



(32) 



and similarly for A^, see Eq. (1291) . Fig. [3] shows the relative variance for the 
propagating state <p as compared to the localized state if) (which was starting 
from the diagonal stationary state). The relative variance undergoes large 
fluctuations, as the variances as such do as well. But in the average, we see 
the expected result. The A(p remains larger than Aip, by about 15 % on time 
average. 



6.2 Ionization properties 



As a further observable, we consider the degree of ionization which, as stated 
above, is a sensitive quantity to probe the effect of SIC. Remind that absorb- 
ing bounds are applied for a proper handling of ionization. We use again an 
instantaneous boost to simulate a very short laser pulse. Fig. |4] shows the re- 
sults for an atom of 3 electrons, within the ID model and for two different 
initial boosts, comparing the TDHF benchmark with TDLDA and TDSIC. 
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Boost=0.2 



Boost=0.3 




5 10 15 20 25 5 10 15 20 25 



Time (fs) Time (fs) 

Fig. 4. Number of escaped electrons, as function of time, for an ID atom of 3 
electrons, initially excited by two different values of boost, as indicated. Three 
time-dependent schemes are compared : TDLDA, TDHF and TDSIC. The model 
parameters are a = 0.5 ao, b = 0.5 ao, R = 0, and z = 0.5. 

The deposited energy with a boost of 0.2 / 0.3 Ha represents 60 % / 134 % of 
the LDA ionization energy (0.100 Ha) and 21 % / 47 % of the SIC ionization 
energy (0.133 Ha). Therefore we expect a strong ionization overestimation for 
LDA. It is obvious that TDSIC comes much closer to the benchmark (TDHF) 
than TDLDA. We checked various other ID molecular systems and found 
similar results confirming that TDSIC recovers nicely the proper ionization 
features. 



6. 3 Results from 3D calculations 

Finally, we want to check the effect of SIC in a realistic 3D situation. To that 
end, we consider the Nas cluster which has a non-symmetric planar structure 
(see insert in Fig. [6]) and whose electron cloud is weakly bound and so pro- 
vides a critical test case for formal developments This cluster contains 
altogether five valence electrons which are active in the low frequency irradia- 
tion processes. The core electrons of the Na atoms are much more bound and 
are eliminated by using pseudo-potentials. For the electronic exchange and 
correlations, we employ the energy-density functional of [39]. The laser exci- 
tation is again simulated by an instantaneous boost. For the further details of 
the 3D calculations see, e.g., [9T3"5] . The first principle result to be mentioned 
is that the newly developed solution scheme for TDSIC runs smoothly also for 
the full 3D case. 

Fig. [5] compares the time evolution of ionization between TDLDA and TDSIC, 
for an initial boost of 0.125 Ha, whose deposited energy represents 149 % of 
the LDA ionization energy (0.105 Ha) and 90 % of the SIC ionization energy 
(0.172 Ha). We see a similarly dramatic effect from SIC on the ionization. 
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Time (fs) 

Fig. 5. Time evolution of the number electrons escaped from the Nas cluster, com- 
puted in full 3D, The system was initially excited by an instantaneous boost. TDSIC 
and TDLDA are compared. 



TDSIC produces less in accordance with the fact that SIC enhances the IP 
from the LDA value of 0.105 Ha to the SIC value 0.172 Ha. 
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Fig. 6. Nas initially excited by an instantaneous boost and computed in full 3D, 
where TDSIC and TDLDA are compared. Left : Time evolution of the dipole spectra 
in x direction. Right : Optical absorption spectrum (in arbitrary units) derived 
thereof [32] ; the insert shows the geometry of Nas • 

Fig. [6] compares the time evolution of the dipole oscillation as such (left) and 
the subsequent optical response (right). The time evolution (left panel) shows 
that the initial oscillations are very much the same for TDLDA and TDSIC. 
Some deviations develop in the further course of propagation which emerges 
from different interferences with particle-hole states. This lets us expect that 
the Mie plasmon position is not affected by TDSIC, but that detailed fragmen- 
tation pattern may be different. These two features are indeed nicely found in 
the optical absorption spectra (right panel). Both results, reduced ionization 
together with little influence on the dominant Mie plasmon excitation was 
also found in earlier studies using a simplified SIC scheme, time- dependent 
average-density SIC [21133] . 
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While Nas is a soft system, as an alternative example of a strongly bound 
system, we consider the case of a C atom. The IP is now 0.224 Ha for LDA, 
whereas it is enhanced to 0.452 Ha for SIC. Then TDSIC produces much 
less ionization than TDLDA when the C atom is excited by an instantaneous 
boost, similarly to the case of Na5 (see Fig. [5]). We also compute the optical 
response spectrum for C, although this observable is less relevant for this non- 
metallic example than for Nas. We observe (not shown here) a shift between 
the SIC peak and the LDA peak. SIC also makes an effect on optical absorption 
because the deeper binding restricts the dipole oscillations more tightly than 
in case of LDA. 



7 Conclusion 

In that paper, we have investigated the time-dependent self-interaction correc- 
tion method (TDSIC) which augments the time- dependent local-density ap- 
proximation (TDLDA) by a self-interaction correction (SIC). That correction 
becomes crucial when aiming at the description of highly dynamical processes 
where electron emission plays a role. But SIC raises problems because the 
resulting one-body Hamiltonian becomes explicitly orbital dependent which, 
in turn, can destroy the necessary orthonormality of the single-particle states. 
One needs to add an explicit constraint on orthonormality. This leads to a 
"symmetry condition", in other words, to the (plausible) condition that the 
SIC mean-field Hamiltonian is hermitian within the space of occupied states. 
An implementation of that involved condition in TDSIC has been achieved by 
dealing with two different sets of single-particle orbitals : Propagating orbitals 
which are carried forth by standard mean-field stepping methods and local- 
izing (or SIC) orbitals which are used to compute the SIC mean-field. The 
relation between the two sets is established by a unitary transformation which 
is determined such that the localizing set satisfies the crucial symmetry condi- 
tion at each time. The newly developed representation of TDSIC constitutes 
a formally consistent and numerically reliable scheme. Crucial conservation 
laws (energy, orthonormality, zero-force theorem) are all obeyed formally and 
in practical calculations. 

First tests have been performed in a one-dimensional model for a molecule. 
By the standard rules of LDA, we have developed a LDA functional for ex- 
change only. That allowed direct comparison with (time-dependent) Hartree- 
Fock (TDHF) as a benchmark. The dynamical evolution was initiated by ex- 
citing the ground state with a very short laser pulse. The pulse was idealized 
as an instantaneous boost and we computed for various excitation strengths 
to check the robustness of the results. It was found that TDSIC compares 
very well with TDHF, while TDLDA overestimates ionization by 50-100%. 
Tests had also been done for fully three dimensional calculations considering 
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a Na5 molecule as test cases. We found again a substantial overestimation of 
ionization for TDLDA, in that case by about 30%. 

Although full 3D calculations have proven to be feasible and stable, the TDSIC 
scheme is numerically costly. We consider IS cLS Sb starting point for further 
developments towards more efficient schemes. The most costly detail is the 
fulfillment of the symmetry condition which basically provides more localized 
single-particle states. This suggests to replace the symmetry condition by a 
direct localization condition which will be numerically less costly. A promising 
option to find it is the time-dependent optimized effective potentials formalism 
[2"Uf2"o] . It has already been used in the static case, leading to what we called a 
"Generalized Slater" potential |44| . The extension to the time dependent case 
is on the way. 
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A Hermiticity of orthonormalization constraints 

The constraint on orthonormality of the single-particle wavefunctions, or cor- 
respondingly unitarity of transformation coefficients, introduces a matrix of 
Lagrangian multipliers and that matrix ought to be hermitian. We prove that 
here for the case of the time- dependent variational principle ( ITT)) . Let us split 
the matrix of Lagrange parameters A into hermitian part \x and anti-hermitian 
part k as A 7j g = yU 7j # + k 7/3 . The variational principle thus becomes = SS for 
the action 




The action S subsequently splits into real and imaginary part where the latter 
reads simply 

%{S}= fdt'^M^ • (A.2) 

Both parts are to be varied independently. Variation of the imaginary part 
yields 

= ^*S{5'} ^^ 7 (r)K 7/3 = V/3,r k 7/3 = . 

7 

This means that we always have ^jS*} = and we deal with a purely real 
action 

The same reasoning applies to all other form of action used for TDSIC paper 
and to the energy functional used for stationary SIC. 



B Alternative derivations of TDSIC 

In this appendix, we will present two alternative derivations of TDSIC. It is 
gratifying to see that different derivations all lead to the same result. The 
alternative routes also shed some new light on the intrinsic properties of TD- 
SIC. 
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B.l The Goedecker method 



First, we will use the Goedecker method of variation [IT]. One starts from 
the SIC energy Esic as given in Eq. (J4]) and defined in term of the the or- 
thonormal set of single particle wavefunctions {ip a }- Thanks to the Lowdin 
orthonormalization method one could equally well expand the {ip a } into 
a set of non-orthogonal functions {i/j a } as 



This is actually another way to constrain the orthonormality of the {ipa}- If 
one assumes that the {ijj a } are not too far from orthonormality, one has to 
first order 



where a a p is a small quantity. One now applies the principle of stationary 
action when varying with respect to the ^* (the non-orthogonal functions). 
This yields 



= 5^ f dt'(£(^|i^) - £sic) • (B.3) 

/3 

Using the chain rule for functional derivatives, we obtain 

The variations with respect to the orthonormal set Sip/3 are similar to those 
used before in the derivations of TDSIC. Its evaluation yields 



The crucial step is now to evaluate the Sip a derivatives. Using Eq. (IB. 21) and 
taking into account that #0* appears not only explicitly in the expansion, but 
also implicitly in the coefficients a al 3, yields 
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= 5 3 (r-r') [5 aP - a aP ] -W ^(r')^(r)5 Q/3 , 

<^*(r') 2^ aK ;m ; 

Note that we can let o~ a p — > after variation. Thus we obtain finally the 
TDSIC equation 

o = cm - Dm - E \^)(^\m. - ^±hi^ a ) (b.4) 

The symmetry condition is recovered by projecting Eq. (1B.4I) on (^g|. 



i?.^ Unitary variation 



Section IB. II uses a variation where orthonormality is explicitly obeyed which 
allows to work without Lagrangian parameters. There is an interesting alter- 
native for such a technique. One deals with a unitary variation according to 
Thouless theorem [16]. Any variation from one Slater state |$) to another 
Slater state |$') can be expressed as 

|$') = exp(ii)|$) , i f = i . (B.5a) 

A variation is a small change and thus the varied states can be obtained from 
linearization. This means 

\5<t>) = L4|$) , = -i($|i t = -i($|i , (B.5b) 

and subsequently in terms of single-particle wavefunctions 

oo 

\5^ a ) = i E \S^ n )A na = i E \5i) P )Ap a + iA ± \S^ a ) , (B.5c) 

n=l /3Gocc 
oo 

(Sip a \ = -i E A«n(^n| = "i E AtfO^I - i(<%*|Ai , (B.5d) 

n=l /3eocc 

i ± = ifl ± + tl ± A , (B.5e) 

where one should be aware of the different ranges of summations, the n running 
over all single-particle space and the (3 only of occupied states. The operator 
A± is the part of the operator leading into space orthogonal to the occupied 
states. Variation then corresponds to a free variation of the matrix elements 
of the hermitian operator A (as long as hermiticity is obeyed). However, the 
limitation to hermiticity means that \Sijj a ) and (5ijj Q \ cannot be varied inde- 
pendently anymore. 
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Applying that variation to the principle of stationary action yields 



o= E ~'^M {m - 1^) - KM dt +hp) a\^p) 

Y. E I -iApnfynl ( ih ~ d t - hp) \ipp) - i(ipp\ (ih dt +hp) \^n)A nll 



/36occ n=l 



J2 E |-i^/3n(^n| (ihd t - hp) \i)p) - i(i/;p\ (ih dt +hp^j \^n)A n p 



PGocc n_L occ 



-i E ^^(^1 (i^ - fy,) |^) - i £ (^| fi^ 9t + M |^ a )4 

/3Gocc /3Gocc 



a/3 



= E E -LA/Sn^nl (i^* ~ V) 1^/3 ) - i(^/s| f ^ <9t +hp) \^ n )A n p \ 
/3Soccn±occ ^ \ / J 

+i E ^/3a(^a|^/3 - ■ 
/3Gocc 

Now, the matrix elements Ag n and v4 n ^ can be varied independently because 
A n p = A*p n e C. Similarly, all elements Ap a can be considered as being inde- 
pendent. This yields the TDSIC equations as 



fl±ihd t ip a = ^-±h a ip a , (B.6a) 
= K a p . (B.6b) 

That, again, reproduces the TDSIC equations (after proper rewriting) together 
with the symmetry condition. 



C LDA exchange energy of the ID model 



One has to evaluate the exchange-correlation energy as a functional of the 
local one-body density : 

Exc[p] = \j d 3 rdV(r(r,r') -p(r)p(r')) (r^r') (C.l) 

For fermionic systems, in the general case : 

r(r,r') = (^|f(r,r'M (C.2) 

where T(r, r') = Z)i>j{^( r ~ r i)^( r ' ~ r j) + <K r ~ r j)$( r ' ~ r i)} is the local two- 
body density matrix (which Pauli effects are included in). In the following, we 
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will focus on the exchange energy only, so that we can limit the tjj to Slater 
determinants. 



In the ID model, one now has to compute : 



T(x,x') = - Y^{ij\T{x,x')\ij), 
hi 



(C.3) 



with : 



f (x, x') = — x i)3(. x ' ~ x j) + $( x ~ Xj)5(x' — Xi)} (C.4) 



i>j 



For a ID free gas, J2i becomes / dk^ (where L is the length of the box) and 

2 f 1 * } 

(ij\T(x, x')\ij) = — I 1 cos [(ki — kj)(x — x')]>. Inserting this in Eq. (1C.3I) . 

one gets : 



T(x, x') 



2^) I dhdkj 



I _ f e Kki~kj)(x-x') 

7 V 



(C.5) 



The first part of the integral gives kp 2 ] the other part is proportional to : 



„ikp(x— x') i p ~ikp(x—x') i 

dk ihix-x') I dk ikj (x-x') = e 



%\x — X' 



-II x — X' 



Collecting these results in ( 1C.5I) . one finds : 



1 — cos [kp{x 

(x — x') 
(C.6) 



T(x, x') = p 2 I 1 - —D [k F (x - x') 



7 



(C.7) 



with the ID free gas density, p = -fjf , and D(x) 



1— COS X 



We now insert (1C.7|) in the ID version of (IC.ip . If the ID potential is 



smoothed Coulomb one, such as 
forward change of variable gives : 



yj (x—x') 2 +a 



(in Hartree units), a straight- 



p 2 f+°° D(k F x) 
pE X (p) = — - / dx- 



7 



7 J-oo \Jx 2 + a (27r) s 



-oo 1 - C0S(^X) 

dx- 



oo X 



l \/x 2 + 



The LDA exchange potential then reads : 



5 



U& K [Q] = jj(pe x {p))\ 



\p=e 



1 f +oo sin^x) 
dx- 



2n J-oo x\J x 2 + a 



(C.9) 
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D Proof of the conservation laws 

D.l Energy conservation 

The proof of energy conservation is straightforward and does not need any 
comments : 

i i 

i i 

D.2 Orthonormality conservation 

The orthonormality should be conserved during time propagation because we 
imposed it in the variation of the action. Nevertheless we will check it ex- 
plicitly Using the TDSIC resulting propagation scheme (13T1) . we see that a 
sufficient condition for the orthonormality to be conserved is that the propa- 
gator is unitary within the space of occupied states, i.e., 



(^|e- Iftsic |^-) = (^-|e- SI %i)* • (D.l) 

This last equation is obviously verified in the occupied subspace (only), be- 
cause of the "weak" hermiticity of hsic in this subspace as given by the sym- 
metry condition (llOcj) . 



D.3 Zero-Force Theorem 



The Zero- Force Theorem (ZFT) states that the kinetic energy plus the electron- 
electron interaction part in the Kohn-Sham mean-field do not change the total 
momentum of the electron cloud which is due to translational symmetry of the 
electron-electron interaction and of the kinetic energy. We formulate that sym- 
bolically as 9f km ' (p) = 0. The proof starts from the Kohn-Sham equations 
where all external fields are dropped 

^ + f> (cl) ll^)=i^ (kin ' el) |^) , (D.2) 
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and where is the Kohn-Sham mean-field part stemming from the elec- 
trons. The time change of total momentum reads 



5 (km,ei) J2(<Pi\p\<Pi) =J2 [(3tV|p|Vi) + (Vi|p|3t kul,eI Vi 

=E [(^t km,d Vi|pvi) + (p^l# m ' el Vi) 

i 

= 4 E f-(tf ( *VlP¥>i) + ( P ^i^ (cl Vi) 

The general form of the ZFT is thus, in {r} representation 



(D.3) 



Ed 



'V^(r) + (r|^ e ^)V^*(r) 



(D.4) 



Now we check whether the TDSIC equation fulfills this theorem, where (r\U^ \ipi 
S« u ia^aUhD A^ill^al 2 }- After simple manipulations on (ID.4[) . one obtains : 

E / d 3 r [(r|C/( el )|^)*V^(r) + (r^l^V^M 

i 

= E / d 3 rf/ LDA ,ei[|^| 2 ]V|^| 2 

= J2[d 3 rU a V Pa (D.5) 
in compact notations. However, since U a is obtained variationally, we have : 



J d 3 r U a Vp a = J d 3 rd 3 



j SE 5p a (r) 
5 p a {r) Sr' 



d d r 



5r' 







(D.6) 



As a consequence, TDSIC does verify the ZFT relation f]D.4[) . as any varia- 
tional scheme should. 
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